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Abstract 

We introduce a numerical method for reconstructing a multidimensional surface using the gradient of the 
surface measured at some values of the coordinates. The method consists of defining a multidimensional 
spline function and minimizing the deviation between its derivatives and the measured gradient. Unlike a 
multidimensional integration along some path, the present method results in a continuous, smooth surface, 
furthermore, it also applies to input data that are non-equidistant and not aligned on a rectangular grid. 
Function values, first and second derivatives and integrals are easy to calculate. The proper estimation of 
the statistical and systematical errors is also incorporated in the method. 
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1. Introduction 

Experiments usually result in a discrete set of datapoints Fi measured at a discrete set of coordinates . 
It is in general useful to interpolate this discrete set to obtain a smooth surface S(q) as a function of the 
coordinates. A more complicated situation is when derivatives dF/dqi are measured, and the surface S(q) 
which fits on this grid of gradients the best is sought for. 

The latter case is commonly encountered in statistical physical problems, e.g. in lattice field theory, where 
using standard Monte-Carlo methods one is not able to measure the (logarithm of the) partition function 
log Z itself, but only its derivatives with respect to the parameters of the theory. These parameters play the 
role of the above coordinates q, and the free energy — log Z as a function of the parameters is what one is 
after. A multidimensional integration is usually carried out in such situations, which, however, can only be 
applied when the measurements reside on a rectangular grid. Furthermore, the result will also depend on the 
integration path and is only guaranteed to be smooth along the path itself. 

The smooth interpolation throughout the whole grid can be obtained by e.g. a multidimensional spline 
function. For a detailed discussion of splines and their application see e.g. Ref. Ref. Q or Ref. Methods 
to determine a two dimensional spline surface upon a rectangular grid given the values of the surface at the 
grid points have been known for a long time (see e.g. Ref. Q), along with algorithms for nonrectangular 
input datasets (see e.g. Ref. i |), using once again the values of the surface at given coordinates. Since 
then spline approximation has received quite a lot of attention and the mathematical basis was studied in 
detail. For recent reviews on these approaches see e.g. Refs. [!,[§]. 

Further methods for multivariate approximation have also been developed, e.g. using rational basis func- 
tions, B-splines, tensor product splines, Powell-Sabin splines, triangulations, genetic algorithms or modified 
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spline techniques. Moreover, based on such methods, various algorithms for spline fitting are also present 



in the literature, see Refs. [lCH26j. Techniques for handling discontinuities (e.g. Ref. (27J) and imposing- 
constraints on the fitted function can also be found (e.g. Refs. 2cl 29|). 



In this paper we present a method which determines the smooth surface S'(q) using the gradient of the 
surface measured at some scattered values of the coordinates, and thus corresponds to a multidimensional 
integration scheme. The algorithm includes an introduction of a set of nodepoints, upon which the multi- 
dimensional spline is determined. This determination is linear and thus straightforward to compute. By a 
systematical variation of the number and position of the nodepoints the method is adapted to the particular 
problem, which is highly important for interpolation in more than one dimensions, as pointed out in Ref. 
Other interpolation schemes for which the parametrization is adapted to the data (like the rational basis 
function approach) are also well suited for the present problem and may be applied in this scope. 

The paper is structured as follows. In section[5]we remind the reader how a spline function in an arbitrary 
number of dimensions D is defined and then in section [3] we show how the fit to the measurements is carried 
out. Since in practice D > 2 is seldom necessary, in order not to complicate the notation we present the 
method in detail for the case of two dimensions. Nevertheless, D > 2 is also straightforward to implement. 
Then we present the algorithm for the systematical placement of the nodepoints, and finally show several 
results on mock data. 

2. Spline definition 

In two dimensions a cubic spline is defined upon a grid {xk,yi} with < k < K, < I < L. The spline 
surface is unambiguously determined by the values that it takes at the nodepoints fk,i — S(xk,yi) (and the 
boundary conditions, which we specify to be "natural", see later). A grid square [xk,Xk+i] x [yi,yi+i] will be 
shortly referred to as {fc, I}. The spline function itself is compactly written a£| 
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The linear equations that determine the spline parameters fk,i from the coefficients for the x and y directions 
(summarized in matrix form as X and Y, respectively) are independent from each other. Inverting these 
equations one obtains for the coefficients 

C V - E E (Y-Xf+M^UJn^ (4) 
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Here the matrix X (Y) only depends on the number K (L) of grid points in the x (y) direction and the widths 
w k i w l )■ T ne c °de for generating the matrix X, together with a sample output is shown in |Appendix A| 
Note that A is a matrix of size A(K — 1) x 4(K — I), but in order to obtain the coefficients we only need the 
first K rows of its inverse - or, even better, the effect of these rows on the vector /. 

Similarly, in arbitrary dimensions D the different directions decouple and thus @ is easily generalized 
to the case of more dimensions. In the following we introduce the fitting method in two dimensions; the 
generalization for higher dimensions is also straightforward to implement. 



^ote that in this formulation the cubic spline contains terms like x 3 y 3 , contrary to other definitions where a bicubic spline 
only has terms x'y J with < i + j < 3. 
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We note here that a compact representation of splines can also be achieved with the help of B-splines 
(for an overview see Refs. [30l432|). In this approach it is not necessary to prescribe particular boundary 
conditions, as for the spline function above (we use "natural" boundary conditions in the definition of X). 



3. Spline fitting 



For any value of the parameters /^j the coefficients C^j - i.e. the whole spline surface S(x, y) - are 
known unambiguously. This way we consistently parameterized the spline surface. Now we want to set the 
spline parameters fk,i such, that derivatives of the spline surface are as close to some previously measured 
values as possible. Note that this way the function S(x, y) will be undetermined upto an overall constant, 
since the translation S(x,y) —¥ S(x,y) + A does not influence the derivatives. This symmetry will be taken 
into account in the following. 

Let us consider N number of points and say that at each point (q^ , Qm ) we measured the derivativ^l in 

the x and y directions: D^} and with errors AD$ and AD$ (m = . . . N - 1). Being "close" can be 
quantified by minimizing 



N-l 



x 2 (f k ,i) = £ 




(5) 



Note that the above % 2 represents a situation where the input data for the derivatives and are 
uncorrelated. However, this is usually not the case and one has to take into account the correlation between 
the measurements. Generalization of the method to include this correlation is presented in |Appcndix B| 

Since S and thus dS/dx and dS/dy are linear in fk,i, this function has a quadratic dependence on the 
parameters fk,i- This enables us to search for the minimum of X 2 (fk,i) 



dx 2 



by solving a system of linear equations 



(6) 



M%l n J k ,l=V ni , n2 , m=0...K-l, n 2 =0...L-l (7) 

Due to the above mentioned translational invariance of the solution, this system of equations is under- 
determined and thus the inverse of M does not exist. This can also be seen by checking that M has a zero 
eigenvalue, corresponding to the eigenvector (1, 1,1,.. .), or, in other words, each row of M adds up to zero. 
Physically this means that one can set e.g. the first element of / to zero, i.e. leave the first column of M. To 
obtain an invertiblc matrix one now has to drop one of its rows, for example the firs|| This way one arrives 
at a matrix A/' of size (KL — 1) x (KL — 1). In the same manner we define V to be the vector composed 
from the last KL — 1 elements of V. One can then complement the solution (M')~ l V with a zero in the first 
element to obtain the final result which satisfies /o,o = 0. 

We remark here that in some cases it can happen that the matrix M is not invertible. An obvious 
example is when there exists a grid square {k,l} in which no measurements reside. The spline function is 
then ambiguous in this grid square and thus M necessarily has a rank smaller than KL — 1. Furthermore, an 
unlucky choice of the coordinates of the measurements q( x ',q( v ' may further reduce the rank of the matrix 
M, but this is very unlikely. By a systematic replacement of the nodepoints most of such situations may be 
filtered out, see section HJ 



2 By measurement we refer to e.g. the lattice determination of the derivatives of log Z at a given value of some bare parameters. 
3 It can be checked that after eliminating the first column, any row of the matrix M can be reproduced by a linear combination 
of the other rows, i.e. it is indeed correct to drop an arbitrary row. 
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In order to explicitly write the elements of the matrix M and the vector V let us analyze the dependence 
of x 2 on the parameters To this end we define k(m) and l(m) as the indices of the grid square that 

contains the mth measurement, i.e. 

(q£\qM)€{k(m),l(m)} (8) 



and define £ m and rj m as the value of the dimensionless coordinate on that grid square corresponding to q 
and q^} (just as in ©): 
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Also, in order to be able to express dS/dx and dS/dy let us define the following matrices: 
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With the matrices 15 and F the expression for x 2 in © can be rewritten as: 
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which implies that in the system of linear equations ([7J to be solved appear 
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Using the actual form of M and V the system of linear equations in ([7]) can be solved^ for With the 
spline parameters the spline coefficients arc also determined through (|4]). 

Since the number of measurements is 2N and we have K-L—l independent spline parameters, the degrees 
of freedom of this fit is given by the difference dof = 2N — K ■ L + 1. This way the freedom corresponding 
to the translational symmetry is transformed out: since /o,o = 0, the result S(x, y) is now set such that 
S(xo,yo) = holds. Note that the spline function can also be transformed easily to satisfy some a priori 
known reference equation S(x r ,y r ) — S r by a simple translation 5 — >• S + (S r — S(x r ,y r )). 

We note that a similar approach to fit a spline function to given gradient information was presented in 
Ref. (H|. 



4. Stable solutions 

The method described in the previous section is bound to give the function S(x, y) for which the sum of 
deviations \ 2 is the smallest. The solution on the other hand also depends on the number and the position 
of the nodepoints, and these have to be tuned appropriately in order to determine the surface that fits best. 

If the number of nodepoints K ■ L is small compared to the number of measurements N, then the reduced 
sum of deviations will be large (x 2 /dof 3> 1) and the spline function S(x, y) may not be a good approximation 



The system of linear equations can be solved using e.g. the Lapack library. 
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to the surface sought for. On the other hand if K ■ L is larg^l, the best fit will become an oscillatory function 
which has the correct derivatives everywhere (i.e. x 2 /dof ~ 1), but is probably not the "right" solution, 
especially when one knows that S should be monotonic|f| in e.g. x. This unwanted feature is a characteristic 
of spline functions even in one dimension. Note that this problem can also be solved using B-splines, where 
the monotonicity of the surface can be guaranteed by imposing linear inequality constraints. 




Figure 1: One dimensional slice of a two dimensional fit. Shown are the input data D^ x ' (red points) and the derivative dS/dx 
of the spline surface obtained using two nodepoint sets (blue and orange bands). The two sets differ only in one gridpoint, 
see Q . On the left side the number of nodepoints in the x direction is K = 10; a slight change in the nodepoints results in no 
visible change in the solution (orange and blue curves are on top of each other). On the other hand, on the right side K = 14, 
and a similar change dramatically distorts the solution. The corresponding values for T> are 0.004 (left) and ~ 10 7 (right). The 
width of the band represents the statistical error. 



Thus we need a measure of how "right" the solution is. A useful way to define this property is to investigate 
how much / changes as the nodepoints are modified, since the oscillatory solutions are unstable even under 
a small change in the nodepoints. This way we can filter out the stable, realistic solutions (see illustration 
on figure [T]). 

We define modified nodepoint sets x^ a ' as 



(a) 
» 



Xk + s, if k = a 
Xk , otherwise 



(14) 



with £ some small number, e.g. e = (xk—1 — Xo)/K/10, and in the same manner for y^>. Now we carry 
out the fitting procedure for each of the modified nodepoint sets, resulting in the modified solutions /w and 
/W. The sum of the relative differences between the original solution / and the modified solutions 
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(15) 



will indeed serve as an indicator of the stability of the fit. If this relative change T> is under a few percents 
then the fit can be considered stable. 



5 Obviously the inequality K ■ L — 1 < 27V should hold otherwise the problem is undcrdetcrmincd. 
6 For example this is the case for the pressure log Z as a function of the temperature. 
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5. Systematics of the method 



The statistical error a stat of the result from the spline fitting method described above can be determined 
using the standard jackknife algorithm, i.e. the system of linear equations ([7]) needs to be solved for each 
jackknife sampl^ll The systematical error on the other hand can be determined by varying the nodepoints 
{xk,Vi}- Based on experience the number of nodepoints may range from M/2 to M; usually an equidistant 
nodepoint-set can already produce a small value for x 2 /dof, but increasing the density of nodepoints in areas 
where the function changes rapidly (i.e. where the measured derivatives are large) can further help to improve 
the fit quality. 

Accordingly, a straightforward way to determine the systematical error is to generate various nodepoint 
sets with different number (and position) of gridpoints. Then, for each set r of the nodepoints the fit is 
carried out resulting in a spline function S T and an indicator G T = (x^ orr /doi) of the fit quality. Results 
which are in the above detailed sense not stable should be filtered out at this point. Then at each point the 
systematical error a sys is determined by 



<7 sys (x,y) = \J (S T (x,y) 2 ) G - (S T (x,y)) 



with 



T T 

Thus the total error can be estimated to be 

(?tot '- 



a sys + a stat 



(16) 
(17) 

(18) 



6. Testing the method against mock data 

We tested the spline fitting method in two dimensions against three sets of mock data. Input to the 
method are the coordinates qff and q$ and the derivatives Dm , with m = . . . N — 1 together with 
10-10 jackknife samples at each m. The derivatives were generated using an original function F{x,y) and 
were scattered for the jackknife samples to have normal distribution with a relative width of A. In table[T] we 
tabulate information about the data: the function F, the number of measurements, the relative error and the 
type of the input distribution (aligned on a rectangular grid or randomly distributed) . The original function 
was chosen such that it resembles the example mentioned earlier: log Z as a function of the temperature and 
some other parameter near a crossover transition (here the role of the temperature is played by the variable 
x). 



fit 


F(x,y) 


N 


A 


input 


1. 
2. 
3. 


(y + 10) 
(V + 2y + 
(2.6y 2 + 2.9y 


■ (2 + tanh(4(.T-4)))(2a; + 3)) 
3) • (1.5 + tanh(4(.T - 4))) (6a; + 
+ 5) • (4 + tanh(3(x - 5))) (3a; - 


3) 

f2) 


400 
1600 
400 


2% 
7% 
2% 


rectangular 
rectangular 
random 



Table 1: Summary of the mock examples. 



After the solution was determined, the fitted surface S(x,y) was compared to the original function and 
to indicate the agreement a weighted sum of deviations 



Is 



Gtotylm ,Qm ) 



N-l 



m— 



(19) 



'Note that this can be performed in a single Lapack call that solves the equations for different vectors on the right hand side. 



() 



was calculated. 

In tabled] we show information about the fits: the minimal value for x 2 /dof and the value of the above 
constructed indicator /3 S . In order to test the method and have a comparison, we also carried out a usual 
two dimensional integration for the rectangular inputs. This was done by integrating a one dimensional 

(x) 

spline of the x-derivatives along the horizontal x-direction upto qin , then the y-derivatives along the vertical 
y-direction upto qmK Then we repeated this procedure in the opposite order, and the difference in the results 
was used to define the systematical error for this method. This way in the same manner as in (|19|) we also 
defined (for rectangular input data) the indicator f3j, which is also shown in the table. In order to compare 
the two methods it is also instructive to study the relative statistical S sta t and relative systematical error S sys 
of the results. In the table we show the average of these quantities S stat and 5 sys for both procedures. 
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X 2 / dof rmn 
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°sys 
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°stat 


°sys 


Pi 


1. 


1.19 


0.14% 


0.27% 


0.47 


0.52% 


0.82% 


0.64 


2. 


1.07 


0.37% 


0.09% 


0.74 


1.66% 


1.36% 


0.35 


3. 


1.33 


0.25% 


0.44% 


0.41 









Table 2: Fit information for the mock examples. 



In all cases results obtained by both methods are consistent with each other within errors. A small value 
for /3 also indicates that the results are indeed good approximations of the original function F(x,y). To 
confirm this we also show in the left side of figure [5] a histogram of the deviation and for fit #2, which 
are well approximated by normal distributions of width < 1. 
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Figure 2: Comparison of the results for fit #2. Histogram of the relative deviation fig and fij between the fitted surface and the 
original function (left side) and histogram of the relative error <5 = J 5^ tat + 5^ ys of the result (right side) . We show the results 
using the spline method with orange, and results obtained by the usual integration with blue. Note that while both procedures 
can be considered as fits with good quality (since the distribution of the deviation fi is narrow) , the relative error is much smaller 
for the spline integration scheme. 

One should however note that the multidimensional spline fitting algorithm results in much smaller statis- 
tical errors than the usual integral method. This is of course due to the fact that in the former case all of the 
measurements (~ N) are taken into account, while for the latter only those lying on the horizontal-vertical 
integration path (~ 2y/~N). The relative error should thus be \/N/V2 times smaller for the spline fitting 
method. The fact that this is indeed realized (see table [2]) shows that the method effectively processes the 
input data. It also is useful to study the distribution of the errors. In the right side of figure[5]a histogram for 
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the total relative error S = yj Sg tat + S] ys is plotted for the case of fit #2. This indicates that for the multi- 
dimensional spline method errors can be an order of magnitude smaller as compared to the naive integration 
procedure. 

We remark that for fit #1 a smaller relative error A in the input results in smaller statistical errors, while 
a large number of input points turns up as a small value for the systematical error. This is also expected, 
since the statistical error is governed by the difference between each jackknife sample, while the systematical 
error depends more on the input data density. We also observe that for the case where the input data were 
randomly distributed, the systematical error dominates over the statistical one. For illustration we plot the 
results S(x,y) for all three examples on figure [3] 




S(x,y) 

8000 
6000 
4000 




Figure 3: The result of the multidimensional spline fitting method for example #1 (left side up), example #2 (right side up) 
and example #3 (down) The original functions F(x,y) were chosen such that they imitate the behaviour of minus one times the 
free energy density in the vicinity of a crossover transition as a function of the temperature (here corresponds to x) and another 
parameter like a bare mass. 

Besides being able to effectively use the increased statistics, the main advantage of the spline fitting 
procedure is that it better estimates the systematics of the result. While for the two-dimensional integration 
only two paths are considered (the vertical- horizontal and the horizontal- vertical paths), the spline method 
is in some sense equivalent to taking into account all possible integration paths at the same time. Our results 
indicate that the contribution coming from these generalized paths cannot be neglected in order to estimate 
the systematics, and the systematic error of the spline method is significantly reduced as compared to the 
naive two-dimensional integration. 



7. Summary 

In this paper we presented an algorithm that determines a smooth hypersurface 5(q) using the measured 
gradient of the surface. The algorithm can thus be applied as a multidimensional integration method, which 
is often useful if one is interested in a continuous approximation of a function of more than one variables. 
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We demonstrated the method on several examples that resemble a situation of the pressure near an analytic, 
crossover-like phase transition, as a function of two parameters. The method, however can easily be imple- 
mented for an arbitrary number of dimensions. Our examples show that usual integration of the gradients 
along a single path is ineffective and does not give a good estimate on the systematics. The proposed method 
on the other hand is capable of taking into account every possible integration path and thus has the ad- 
vantage that it both decreases the statistical error and better estimates the systematical error of the result. 
The method is not restricted to input data that resides on a rectangular grid, but is capable of handling the 
situation when the gradient of the function is measured at scattered values of the parameters q. 

The new method is primarily applicable in statistical physics, e.g. in lattice field theory. Here it constitutes 
an essential improvement as compared to the typically used methods like the conventional integration. This 
scope of application also motivated the example surfaces of section [51 which resemble continuous phase 
transitions in the space of some bare parameters. In particular, the present method was utilized to obtain 



results for the pressure log Z using the lattice determined derivatives of log Z in Ref . [34j . 
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Appendix A. The spline matrix 

In the following we list the code that generates the matrix X: 
/* Input */ 

int K; // Number of grid points in the given direction} 

int n_exp; // Number of exponents; this works for n_exp=4 or 5 
double *W; // Array of size (K-l) containing the grid widths 

/* Output */ 

double *X // Array of size (n_exp* (K-l) ) * (n_exp* (K-l) ) , result 

/* Auxiliary variables */ 
int k , 1 , p ; 
int eqno = ; 

int fac[4][5] = {{1 , 1 , 1 , 1 , 1} , {0 , 1 ,2 , 3,4} , {0 ,0 , 2 ,6, 12} ,{0 ,0 ,0 ,6 ,24}} ; 
int n = n_exp * (K-l) ; 

/* fill X with zeros */ 

for (k=0; k<n; k++) for (1=0; Kn; 1++) { 
X[k*n+1] =0.0; 

} 

/* fix spline value at K points */ 
for (1=0; KK-1; 1++) { 

X[eqno*n + n_exp*l] = 1.0; 

eqno++ ; 

} 

for (k=0; k<n_exp; k++) X[eqno*n + n_exp*(K-2)+k] = 1.0; 
eqno++ ; 

/* fix continuity of zeroth, first, second ... (n_exp-2)-th derivative of spline */ 
for (1=0; KK-2; 1++) for (p=0; p<n_exp-l; p++) { 

for (k=0; k<n_exp; k++) X[eqno*n + n_exp*l + k] = fac[p][k] / pow(W[l],p); 

X[eqno*n + n_exp*(l+l) + p] = -fac[p][p] / pow(W[l+l] ,p) ; 

eqno++ ; 

} 

/* (n_exp-2)-th derivative goes to zero at the ends */ 

X[eqno*n + n_exp-2] = f ac [n_exp-2] [n_exp-2] / pow(W [0] ,n_exp-2) ; 

eqno++; 

for (k=0; k<n_exp; k++) 

X[eqno*n + (K-2)*n_exp + k] = f ac [n_exp-2] [k] / pow(W [K-2] ,n_exp-2) ; 
eqno++; 

/* for n_exp=5 set second derivative to zero at left end */ 
if (n_exp == 5) { 

X[eqno*n + n_exp-3] = f ac [n_exp-3] [n_exp-3] / pow(W [0] ,n_exp-3) ; 

eqno++ ; 

} 
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Using the input values K=5, n_exp=4 and the vector of widths iu, the output from the above code for X is 
the following: 
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Appendix B. Correlated case 



If the measured derivatives D x and D y are correlated, the \ 2 function contains additional terms. Usually 
measurements at different m values are independent, but the correlation of D x (m) and D y (rn) at the same m 
is not negligible, since e.g. in a lattice calculation these are determined using the same configurations. This 
leads to the following xl 



< 2 ■ 

Vcorr 



\ 2 

/Vcorr 



M-l 

E 

m=0 



dS_ 

'Moo V dx ' 
2Q- 1 (- 



D x {m) + 



- Mm)) (f - D y (m)) + Q^ )n - D y (m) 



(B.l) 



where Q( m ) is the 2x2 correlation matrix consisting of the correlators of the two measured derivatives at 
the mth point. The corresponding system of linear equations has the same form as ([7]); only now on the left 
and right hand side enter 



N-l 



lvl n!,n 2 — P(m) 00 ™ "i ^ V (™)oi m m m m ' m m 



m=0 
N-l 



Yn 



E [0(i ^(m)C'" 2 + 2Q~l )n {D x {m)FZr 2 + E^D y {m)) + Q^D^F^ 



m=0 



(B.2) 



(B.3) 
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